################################################
# Analysis for:
#
# Ferguson, Grant, Kathryn Haglin, and Soren Jordan. "The Personality and Politics of Cryptocurrency
#  Investors." At American Politics Research. 
#  Partially funded by a Grant-in-Aid of Research, Artistry and Scholarship for $27,500. 
#  DOI: 10.1177/1532673X231220653
#
# Required files: FHJ-APR-2023.csv
#
################################################

library(tidyverse)
library(broom)
library(nnet)

data <- read_csv("/Users/scj0014/Myfiles/Research/Cryptocurrency/Personality/To website/Replication/FHJ-APR-2023.csv")

names(data)

# educ
#	1 = No high school
#	2 = High school
#	3 = Some college
# 	4 = 2-year degree
#	5 = 4-year degree
#	6 = Post-grad

# age
#	<> = Numeric age

# male (YouGov does not have third categories)
#	0 = Female
#	1 = Male

# black
#	0 = Not Black
#	1 = Black

# hispanic
#	0 = Not Hispanic
#	1 = Hispanic

# asian
#	0 = Not Asian
#	1 = Asian

# white
#	0 = Not white
#	1 = White

# married
#	0 = Not married
#	1 = Married

# religimportance
# 	1 = Not at all important
#	2 = Not too important
#	3 = Somewhat important
#	4 = Very important

# income
# 	1 = Less than $10,000
#	2 = $10,000 - $19,999
#	3 = $20,000 - $29,999
#	4 = $30,000 - $39,999
#	5 = $40,000 - $49,999
#	6 = $50,000 - $59,999
#	7 = $60,000 - $69,999 
#	8 = $70,000 - $79,999
#	9 = $80,000 - $99,999
#	10 = $100,000 - $119,999
#	11 = $120,000 - $149,999
#	12 = $150,000 - $199,999
#	13 = $200,000 - $249,999
#	14 = $250,000 - $349,999
#	15 = $350,000 - $499,999
#	16 = $500,000 or more

# stockowner
#	0 = Does not own stocks
#	1 = Owns stocks

# popdensity
#	1 = Rural area
#	2 = Small town
#	3 = Suburban area
#	4 = Smaller city
#	5 = Big city

# bornagain
#	0 = Is not born-again
#	1 = Is born-again

# catholic
#	0 = Is not Catholic
#	1 = Is Catholic

# openexp. Average of two 1-7 scales (open to new experiences; not conventional)
# 	1 = Strongly disagree (I am not those traits / open to experience)
# 	...
#	7 = Strongly agree (I am those traits / open to experience)

# conscientiousness. Average of two 1-7 scales (dependable; not disorganized)
# 	1 = Strongly disagree (I am not those traits / conscientious)
# 	...
#	7 = Strongly agree (I am those traits / conscientious)

# extraversion. Average of two 1-7 scales (extraverted; not reserved)
# 	1 = Strongly disagree (I am not those traits / extraverted)
# 	...
#	7 = Strongly agree (I am those traits / extraverted)

# agreeableness. Average of two 1-7 scales (sympathetic; not quarrelsome)
# 	1 = Strongly disagree (I am not those traits / agreeable)
# 	...
#	7 = Strongly agree (I am those traits / agreeable)

# emostable. Average of two 1-7 scales (calm; not anxious)
# 	1 = Strongly disagree (I am not those traits / emotionally stable)
# 	...
#	7 = Strongly agree (I am those traits / emotionally stable)

# n2e. Average of two 0-1 scales (opinions versus others; think about things)
# 	0 = Low need to evaluate
# 	...
#	1 = High need to evaluate

# n2e. Average of two 0-1 scales (responsibiltiy for thinking; simple vs. complex problems)
# 	0 = Low need to evaluate
# 	...
#	1 = High need to evaluate

# machiavellianism. Average of four 1-7 scales (manipulate others; lie to get way; use flattery; exploit others)
# 	1 = Strongly disagree (I am not those traits / Machiavellian)
# 	...
#	7 = Strongly agree (I am not those traits / Machiavellian)

# narcissism. Average of four 1-7 scales (others admire me; attention to self; get favors; seek prestige)
# 	1 = Strongly disagree (I am not those traits / narcissistic)
# 	...
#	7 = Strongly agree (I am not those traits / narcissistic)

# psychopathy. Average of four 1-7 scales (lack remorse; don't care about morality; callous; cynical)
# 	1 = Strongly disagree (I am not those traits / psychopathic)
# 	...
#	7 = Strongly agree (I am not those traits / psychopathic)

# spendserve
#	1 = Government should strongly decrease overall spending, even if it means fewer services
#	2 = Government should somewhat decrease overall spending, even if it means fewer services
#	3 = Government spends about the right amount and provides about the right number of services
#	4 = Government should somewhat increase overall services, even if it means spending more
#	5 = Government should strongly increase overall services, even if it means spending more

# personalinflation. Personal hardship due to inflation
#	0 = No personal hardship
#	...
#	10 = Maximum personal hardship

# inflationexpectations
#	0 = Prices will go down / stay the same
#	1 = Prices will go up

# trustgov
#	1 = Never
#	2 = Some of the time
#	3 = About half of the time
#	4 = Most of the time
#	5 = Always

# conspiratorialthinking. Average of four 1-5 scales (plots; few people run the country; who runs not known; big events)
#	1 = Strongly disagree (not conspiratorial)
#	2 = Somewhat disagree
#	3 = Neither agree nor disagree
#	4 = Somewhat agree
#	5 = Strongly agree (conspiratorial)

# moraltrad. Average of four 1-5 scales (family ties; new lifestyles; world changes morality; tolerate other standards)
#	1 = Strongly disagree (not morally traditional)
#	2 = Somewhat disagree
#	3 = Neither agree nor disagree
#	4 = Somewhat agree
#	5 = Strongly agree (morally traditional)

# egal. Average of four 1-5 scales (equal opportunity; worry less about equality; others have more chance; equality less problems)
#	1 = Strongly disagree (not morally egalitarian)
#	2 = Somewhat disagree
#	3 = Neither agree nor disagree
#	4 = Somewhat agree
#	5 = Strongly agree (morally egalitarian)

# libertarianvalues. average of egal and moraltrad (which ranges from 1-5), but inverted. So low on egal/moraltrad; high on this
#	1 = Strongly morally traditional <and> egalitarian: NOT libertarian
#	...
#	5 = Strongly not morally traditional <or> egalitarian: libertarian

# symbideologycons
#	1 = Extremely liberal
#	2 = Liberal
#	3 = Slightly liberal
#	4 = Moderate, middle of the road
#	5 = Slightly conservative
#	6 = Conservative
#	7 = Extremely conservative

# pid7sr
#	1 = Strong Democrat
#	2 = Not very strong Democrat
#	3 = Lean Democrat
#	4 = Independent
#	5 = Lean Republican
#	6 = Not very strong Republican
#	7 = Extremely Republican

# cryptoowner
#	0 = Not a cryptocurrency owner
#	1 = Cryptocurrency owner


# Create 0-1 versions so we can compare coefficient plots to each other
range01 <- function(x) { 
	(x - min(x, na.rm = TRUE)) / (max(x, na.rm = TRUE) - min(x, na.rm = TRUE))
}

data01 <- data.frame(apply(data, 2, range01))



###############################################################################################
# Figure 1. Symbolic ideology of cryptocurrency owners
###############################################################################################

ggplot(data = data[data$cryptoowner == 1,], aes(x = symbideologycons)) +
	geom_bar() +
	labs(y = "Number of Respondents", x = "Self-Placed Ideology") +
	scale_x_continuous(labels = c("Extremely\nLiberal", "Liberal", "Somewhat\nLiberal", "Moderate",
		"Somewhat\nConservative", "Conservative", "Extremely\nConservative"), breaks = 1:7) + 
	theme_bw()
dev.off()


###############################################################################################
# Figure 2. Main models: predictors of cryptocurrency ownership
###############################################################################################

## Base model
base.model <- glm(cryptoowner ~ age + male + black + hispanic + asian + educ + married + popdensity +
								bornagain + catholic + religimportance + 
								income + stockowner + # personalinflation + inflationexpectations +
								# spendserve + trustgov + conspiratorialthinking + symbideologycons + libertarianvalues +
								openexp + conscientiousness + extraversion + agreeableness + emostable + 
								n2e + n4c + narcissism + machiavellianism + psychopathy, family = binomial(link = "logit"), data = data01)

summary(base.model)

# Transform to tidy dataframe, lose the intercept for plotting
base.model.plot <- tidy(base.model)
base.model.plot <- base.model.plot[!(base.model.plot$term %in% "(Intercept)"),]

# Pretty labels for plotting
the.base.predictors <- c("Age", "Male", "Black", "Hispanic", "Asian", "Education", "Married", "Population\nDensity",
							"Born-Again", "Catholic", "Importance of\nReligion", 
							"Income", "Stock Owner", # "Personal Inflation\nHardship", "Inflation\nExpectations",
							# "Pro-Spending\nfor Services", "Trust in\nGovernment", "Conspiratorial\nThinking",
							# "Symbolic\nIdeology", "Libertarian\nValues",
							"Openness to\nExperience", "Conscientiousness", "Extraversion", "Agreeableness", "Emotional\nStability",
							"Need to Evaluate", "Need for Cognition", "Narcissism", "Machiavellianism", "Psychopathy")

# Factored labels in correct order for plot
base.model.plot$Predictor <- factor(the.base.predictors, 
								levels = rev(the.base.predictors)) # ordered in the regression, will reverse on plotting with flip

# Categories of predictors for facet
base.model.plot$Category <- as.factor(c(rep("Demographics", 11), rep("Economics", 2), 
								rep("Politics", 0), rep("Personality", 10)))
								
# Manually assigned, so double check
print(base.model.plot, n = 50)

# Error bars: 0.05 and 0.10 since so few owners
base.model.plot$lower.90 <- base.model.plot$estimate - 1.65*base.model.plot$std.error
base.model.plot$upper.90 <- base.model.plot$estimate + 1.65*base.model.plot$std.error
base.model.plot$lower.95 <- base.model.plot$estimate - 1.96*base.model.plot$std.error
base.model.plot$upper.95 <- base.model.plot$estimate + 1.96*base.model.plot$std.error

# Assign model type for color
base.model.plot$Model <- "Base"


## Full model
full.model <- glm(cryptoowner ~ age + male + black + hispanic + asian + educ + married + popdensity +
								bornagain + catholic + religimportance + 
								income + stockowner + personalinflation + inflationexpectations +
								spendserve + trustgov + conspiratorialthinking + symbideologycons + libertarianvalues +
								openexp + conscientiousness + extraversion + agreeableness + emostable + 
								n2e + n4c + narcissism + machiavellianism + psychopathy, 
								family = binomial(link = "logit"), data = data01)

summary(full.model)

# Transform to tidy dataframe, lose the intercept for plotting
full.model.plot <- tidy(full.model)
full.model.plot <- full.model.plot[!(full.model.plot$term %in% "(Intercept)"),]

# Pretty labels for plotting
the.full.predictors <- c("Age", "Male", "Black", "Hispanic", "Asian", "Education", "Married", "Population\nDensity",
							"Born-Again", "Catholic", "Importance of\nReligion", 
							"Income", "Stock Owner", "Personal Inflation\nHardship", "Inflation\nExpectations",
							"Pro-Spending\nfor Services", "Trust in\nGovernment", "Conspiratorial\nThinking",
							"Symbolic\nIdeology", "Libertarian\nValues",
							"Openness to\nExperience", "Conscientiousness", "Extraversion", "Agreeableness", "Emotional\nStability",
							"Need to Evaluate", "Need for Cognition", "Narcissism", "Machiavellianism", "Psychopathy")

# Factored labels in correct order for plot
full.model.plot$Predictor <- factor(the.full.predictors, 
								levels = rev(the.full.predictors)) # ordered in the regression, will reverse on plotting with flip

# Categories of predictors for facet
full.model.plot$Category <- as.factor(c(rep("Demographics", 11), rep("Economics", 4), 
								rep("Politics", 5), rep("Personality", 10)))
								
# Manually assigned, so double check
print(full.model.plot, n = 50)

# Error bars: 0.05 and 0.10 since so few owners
full.model.plot$lower.90 <- full.model.plot$estimate - 1.65*full.model.plot$std.error
full.model.plot$upper.90 <- full.model.plot$estimate + 1.65*full.model.plot$std.error
full.model.plot$lower.95 <- full.model.plot$estimate - 1.96*full.model.plot$std.error
full.model.plot$upper.95 <- full.model.plot$estimate + 1.96*full.model.plot$std.error

# Assign model type for color
full.model.plot$Model <- "Full"

# Combine
both.model.plot <- bind_rows(base.model.plot, full.model.plot)

ggplot(data = both.model.plot, aes(x = Predictor, y = estimate, color = Model)) +
	geom_hline(yintercept = 0, lwd = 1.25, color = "grey80") + 
	geom_point(size = 3, position = position_dodge(width = 0.5)) +
	scale_color_manual(values = c("#DF65B0", "#67001F"), breaks = c("Full", "Base")) +
	geom_errorbar(aes(ymin = lower.95, ymax = upper.95), lwd = 0.5, width = 0, position = position_dodge(width = 0.5)) +
	geom_errorbar(aes(ymin = lower.90, ymax = upper.90), lwd = 1.1, width = 0, position = position_dodge(width = 0.5)) +
	facet_wrap(~ Category, scales = "free_y") +
	labs(y = "Estimate (All Predictors Scaled from 0-1)") +
	coord_flip() +
	theme_bw()
dev.off()






###############################################################################################
# Figure 3. Predicted probabilities of cryptocurrency ownership from full model
###############################################################################################

# Establish baseline
profile.base <- data.frame(age = as.numeric(names(table(data01$age))[names(table(data$age)) == 50]), # convert to 01
				male = 0,
				black = 0,
				hispanic = 0,
				asian = 0, 
				educ = as.numeric(names(table(data01$educ))[names(table(data$educ)) == 4]), # convert to 01
				married = 0,
				popdensity = as.numeric(names(table(data01$popdensity))[names(table(data$popdensity)) == 3]), # convert to 01
				bornagain = 0,
				catholic = 0,
				religimportance = as.numeric(names(table(data01$religimportance))[names(table(data$religimportance)) == 3]), # convert to 01
				income = as.numeric(names(table(data01$income))[names(table(data$income)) == 9]),  # convert to 01
				stockowner = 1,
				personalinflation = mean(data01$personalinflation, na.rm = TRUE),
				inflationexpectations = mean(data01$inflationexpectations, na.rm = TRUE),
				spendserve = mean(data01$spendserve, na.rm = TRUE),
				trustgov = mean(data01$trustgov, na.rm = TRUE),
				conspiratorialthinking = mean(data01$conspiratorialthinking, na.rm = TRUE),
				symbideologycons = mean(data01$symbideologycons, na.rm = TRUE),
				libertarianvalues = mean(data01$libertarianvalues, na.rm = TRUE),
				openexp = mean(data01$openexp, na.rm = TRUE),
				conscientiousness = mean(data01$conscientiousness, na.rm = TRUE),
				extraversion = mean(data01$extraversion, na.rm = TRUE),				
				agreeableness = mean(data01$agreeableness, na.rm = TRUE),
				emostable = mean(data01$emostable, na.rm = TRUE),
				n2e = mean(data01$n2e, na.rm = TRUE),
				n4c = mean(data01$n4c, na.rm = TRUE),				
				narcissism = mean(data01$narcissism, na.rm = TRUE),				
				machiavellianism = mean(data01$machiavellianism, na.rm = TRUE),				
				psychopathy = mean(data01$psychopathy, na.rm = TRUE)
)

				
profile.base.fit <- predict(full.model, profile.base, type = "response", se.fit = TRUE)
profile.base.data <- data.frame(Predictor = "Baseline\nRespondent", value = profile.base.fit$fit,
				value.lower = profile.base.fit$fit - 1.96*profile.base.fit$se.fit,
				value.upper = profile.base.fit$fit + 1.96*profile.base.fit$se.fit,
				Category = "Baseline")

## Update for each. 
# Shift baseline from 50 to 30 for age
profile.age <- profile.base
profile.age$age <- as.numeric(names(table(data01$age))[names(table(data$age)) == 30])
profile.age.fit <- predict(full.model, profile.age, type = "response", se.fit = TRUE)
profile.age.data <- data.frame(Predictor = "Age: 50 -> 30 years", value = profile.age.fit$fit,
				value.lower = profile.age.fit$fit - 1.96*profile.age.fit$se.fit,
				value.upper = profile.age.fit$fit + 1.96*profile.age.fit$se.fit,
				Category = "Demographics")

# Shift baseline from 0 to 1 for male
profile.male <- profile.base
profile.male$male <- 1
profile.male.fit <- predict(full.model, profile.male, type = "response", se.fit = TRUE)
profile.male.data <- data.frame(Predictor = "Male: Female -> Male", value = profile.male.fit$fit,
				value.lower = profile.male.fit$fit - 1.96*profile.male.fit$se.fit,
				value.upper = profile.male.fit$fit + 1.96*profile.male.fit$se.fit,
				Category = "Demographics")

# Shift baseline from 1 to 0 for stockowner			
profile.nostocks <- profile.base
profile.nostocks$stockowner <- 0
profile.nostocks.fit <- predict(full.model, profile.nostocks, type = "response", se.fit = TRUE)
profile.nostocks.data <- data.frame(Predictor = "Stockowner: Yes -> No", value = profile.nostocks.fit$fit,
				value.lower = profile.nostocks.fit$fit - 1.96*profile.nostocks.fit$se.fit,
				value.upper = profile.nostocks.fit$fit + 1.96*profile.nostocks.fit$se.fit,
				Category = "Economics")

# Shift baseline from 0 to 1 for hispanic
profile.hispanic <- profile.base
profile.hispanic$hispanic <- 1
profile.hispanic.fit <- predict(full.model, profile.hispanic, type = "response", se.fit = TRUE)
profile.hispanic.data <- data.frame(Predictor = "Hispanic: No -> Yes", value = profile.hispanic.fit$fit,
				value.lower = profile.hispanic.fit$fit - 1.96*profile.hispanic.fit$se.fit,
				value.upper = profile.hispanic.fit$fit + 1.96*profile.hispanic.fit$se.fit,
				Category = "Demographics")

# Shift baseline from 0 to 1 for black			
profile.black <- profile.base
profile.black$black <- 1
profile.black.fit <- predict(full.model, profile.black, type = "response", se.fit = TRUE)
profile.black.data <- data.frame(Predictor = "Black: No -> Yes", value = profile.black.fit$fit,
				value.lower = profile.black.fit$fit - 1.96*profile.black.fit$se.fit,
				value.upper = profile.black.fit$fit + 1.96*profile.black.fit$se.fit,
				Category = "Demographics")

# Shift baseline from mean to maximum for openness to experience
profile.openexp <- profile.base
profile.openexp$openexp <- 1   # max on 01
profile.openexp.fit <- predict(full.model, profile.openexp, type = "response", se.fit = TRUE)
profile.openexp.data <- data.frame(Predictor = "Openness to Experience:\nMean -> Maximum", value = profile.openexp.fit$fit,
				value.lower = profile.openexp.fit$fit - 1.96*profile.openexp.fit$se.fit,
				value.upper = profile.openexp.fit$fit + 1.96*profile.openexp.fit$se.fit,
				Category = "Personality")

# Shift baseline from mean to minimum for conscientiousness
profile.conscientiousness <- profile.base
profile.conscientiousness$conscientiousness <- 0   # min on 01
profile.conscientiousness.fit <- predict(full.model, profile.conscientiousness, type = "response", se.fit = TRUE)
profile.conscientiousness.data <- data.frame(Predictor = "Conscientiousness:\nMean -> Minimum", value = profile.conscientiousness.fit$fit,
				value.lower = profile.conscientiousness.fit$fit - 1.96*profile.conscientiousness.fit$se.fit,
				value.upper = profile.conscientiousness.fit$fit + 1.96*profile.conscientiousness.fit$se.fit,
				Category = "Personality")

# Shift baseline from mean to maximum for narcissism
profile.narcissism <- profile.base
profile.narcissism$narcissism <- 1   # max on 01
profile.narcissism.fit <- predict(full.model, profile.narcissism, type = "response", se.fit = TRUE)
profile.narcissism.data <- data.frame(Predictor = "Narcissism: Mean -> Maximum", value = profile.narcissism.fit$fit,
				value.lower = profile.narcissism.fit$fit - 1.96*profile.narcissism.fit$se.fit,
				value.upper = profile.narcissism.fit$fit + 1.96*profile.narcissism.fit$se.fit,
				Category = "Personality")
	
# Shift baseline from mean to minimum for personal inflation hardship			
profile.personalinflation.min <- profile.base
profile.personalinflation.min$personalinflation <- 0   # min on 01
profile.personalinflation.min.fit <- predict(full.model, profile.personalinflation.min, type = "response", se.fit = TRUE)
profile.personalinflation.min.data <- data.frame(Predictor = "Personal Inflation Hardship:\nMean -> Minimum", value = profile.personalinflation.min.fit$fit,
				value.lower = profile.personalinflation.min.fit$fit - 1.96*profile.personalinflation.min.fit$se.fit,
				value.upper = profile.personalinflation.min.fit$fit + 1.96*profile.personalinflation.min.fit$se.fit,
				Category = "Economics")

# Shift baseline from mean to maximum for personal inflation hardship
profile.personalinflation.max <- profile.base
profile.personalinflation.max$personalinflation <- 1   # max on 01
profile.personalinflation.max.fit <- predict(full.model, profile.personalinflation.max, type = "response", se.fit = TRUE)
profile.personalinflation.max.data <- data.frame(Predictor = "Personal Inflation Hardship:\nMean -> Maximum", value = profile.personalinflation.max.fit$fit,
				value.lower = profile.personalinflation.max.fit$fit - 1.96*profile.personalinflation.max.fit$se.fit,
				value.upper = profile.personalinflation.max.fit$fit + 1.96*profile.personalinflation.max.fit$se.fit,
				Category = "Economics")
				
# Shift baseline from mean to minimum for spending versus services
profile.spendserve.min <- profile.base
profile.spendserve.min$spendserve <- 0   # min on 01
profile.spendserve.min.fit <- predict(full.model, profile.spendserve.min, type = "response", se.fit = TRUE)
profile.spendserve.min.data <- data.frame(Predictor = "Pro-Spending for Services:\nMean -> Minimum", value = profile.spendserve.min.fit$fit,
				value.lower = profile.spendserve.min.fit$fit - 1.96*profile.spendserve.min.fit$se.fit,
				value.upper = profile.spendserve.min.fit$fit + 1.96*profile.spendserve.min.fit$se.fit,
				Category = "Politics")

# Shift baseline from mean to maximum for spending versus services
profile.spendserve.max <- profile.base
profile.spendserve.max$spendserve <- 1   # max on 01
profile.spendserve.max.fit <- predict(full.model, profile.spendserve.max, type = "response", se.fit = TRUE)
profile.spendserve.max.data <- data.frame(Predictor = "Pro-Spending for Services:\nMean -> Maximum", value = profile.spendserve.max.fit$fit,
				value.lower = profile.spendserve.max.fit$fit - 1.96*profile.spendserve.max.fit$se.fit,
				value.upper = profile.spendserve.max.fit$fit + 1.96*profile.spendserve.max.fit$se.fit,
				Category = "Politics")

# Shift baseline from mean to maximum for conspiratorial thinking
profile.conspiratorialthinking.max <- profile.base
profile.conspiratorialthinking.max$conspiratorialthinking <- 1   # max on 01
profile.conspiratorialthinking.max.fit <- predict(full.model, profile.conspiratorialthinking.max, type = "response", se.fit = TRUE)
profile.conspiratorialthinking.max.data <- data.frame(Predictor = "Conspiratorial Thinking:\nMean -> Maximum", value = profile.conspiratorialthinking.max.fit$fit,
				value.lower = profile.conspiratorialthinking.max.fit$fit - 1.96*profile.conspiratorialthinking.max.fit$se.fit,
				value.upper = profile.conspiratorialthinking.max.fit$fit + 1.96*profile.conspiratorialthinking.max.fit$se.fit,
				Category = "Politics")

# Combine predicted probabilities together
pred.prob.data <- bind_rows(profile.base.data,
	profile.age.data,
	profile.male.data,
	profile.nostocks.data,
	profile.hispanic.data,
	profile.black.data,
	profile.openexp.data,
	profile.conscientiousness.data,
	profile.narcissism.data,
	profile.personalinflation.min.data,
	profile.personalinflation.max.data,
	profile.spendserve.min.data,
	profile.spendserve.max.data,
	profile.conspiratorialthinking.max.data
)

# Factored labels in correct order for plot
pred.prob.data$Predictor.order <- factor(pred.prob.data$Predictor, levels = rev(c( 
	"Baseline\nRespondent", "Age: 50 -> 30 years", "Male: Female -> Male", "Black: No -> Yes", "Hispanic: No -> Yes",		
	"Personal Inflation Hardship:\nMean -> Minimum", "Personal Inflation Hardship:\nMean -> Maximum", "Stockowner: Yes -> No",
	"Openness to Experience:\nMean -> Maximum", "Conscientiousness:\nMean -> Minimum", "Narcissism: Mean -> Maximum",                   
	"Pro-Spending for Services:\nMean -> Minimum", "Pro-Spending for Services:\nMean -> Maximum", "Conspiratorial Thinking:\nMean -> Maximum")))

ggplot(data = pred.prob.data, aes(x = Predictor.order, y = value)) +
	geom_point(size = 3, color = "#67001F") +
	geom_errorbar(aes(ymin = value.lower, ymax = value.upper), lwd = 0.5, width = 0, color = "#67001F") +
	labs(y = "Predicted Probabiltiy of Cryptocurrency Ownership", x = "Predictor") +
	coord_flip() +
	theme_bw()
dev.off()




###############################################################################################
# Figure 4. Multinomial logit of investment choices
###############################################################################################

# Create four distinct categories: table of ownership of stocks + crypto (add to both versions of dataset)
data$cryptostock <- data01$cryptostock <- ifelse((data$cryptoowner == 0 & data$stockowner == 0), "Neither", 
				ifelse((data$cryptoowner == 0 & data$stockowner == 1), "Just stocks",
					ifelse((data$cryptoowner == 1 & data$stockowner == 0), "Just crypto", "Both")))

table(data$cryptostock, data$cryptoowner)
table(data$cryptostock, data$stockowner)

# Establish ``Neither'' as baseline
data$cryptostock <-  data01$cryptostock <- factor(data$cryptostock, levels = c("Neither", "Just stocks", "Just crypto", "Both"))


table(data01$asian, data01$cryptostock)
# Only to Asian respondents reply ``Just crypto''. The model is not identified for them (separation). 
#  So we must omit this demographic predictor from the model. I am okay with this, though, since it was
#  not a significant predictor in the original regression.

model.multinom <- multinom(cryptostock ~ age + male + black + hispanic + #  asian + 
								educ + married + popdensity +
								bornagain + catholic + religimportance + 
								income + personalinflation + inflationexpectations +
								spendserve + trustgov + conspiratorialthinking + symbideologycons + libertarianvalues +
								openexp + conscientiousness + extraversion + agreeableness + emostable + 
								n2e + n4c + narcissism + machiavellianism + psychopathy, data = data01)

# Grab estimates and standard errors from three separate models for plotting
model.multi.plot <- data.frame(estimate = c(coef(model.multinom)[1,], coef(model.multinom)[2,], coef(model.multinom)[3,]),
						std.error = sqrt(diag(vcov(model.multinom))),
						Model = c(rep("Neither -> Just Stocks", length(coef(model.multinom)[1,])), 
								rep("Neither -> Just Crypto", length(coef(model.multinom)[1,])),
								rep("Neither -> Both", length(coef(model.multinom)[1,]))))

# Pretty labels for plotting
model.multi.plot$Predictor <- model.multinom.predictors <- c("(Intercept)", "Age", 
							"Male", "Black", "Hispanic", # "Asian", 
							"Education", "Married", "Population\nDensity",
							"Born-Again", "Catholic", "Importance of\nReligion", 
							"Income", "Personal Inflation\nHardship", "Inflation\nExpectations",
							"Pro-Spending\nfor Services", "Trust in\nGovernment", "Conspiratorial\nThinking",
							"Symbolic\nIdeology", "Libertarian\nValues",
							"Openness to\nExperience", "Conscientiousness", "Extraversion", "Agreeableness", "Emotional\nStability",
							"Need to Evaluate", "Need for Cognition", "Narcissism", "Machiavellianism", "Psychopathy")

# Factored labels in correct order for plotting
model.multi.plot$Predictor <- factor(model.multi.plot$Predictor, 
								levels = rev(model.multinom.predictors)) # ordered in the regression, will reverse on plotting with flip

model.multi.plot <- model.multi.plot[!(model.multi.plot$Predictor %in% "(Intercept)"),]

# Establish confidence intervals around coefficients
model.multi.plot$lower.90 <- model.multi.plot$estimate - 1.65*model.multi.plot$std.error
model.multi.plot$upper.90 <- model.multi.plot$estimate + 1.65*model.multi.plot$std.error
model.multi.plot$lower.95 <- model.multi.plot$estimate - 1.96*model.multi.plot$std.error
model.multi.plot$upper.95 <- model.multi.plot$estimate + 1.96*model.multi.plot$std.error

# Categories of predictors for facet
model.multi.plot$Category <- as.factor(c(rep("Demographics", 10), rep("Economics", 3), 
								rep("Politics", 5), rep("Personality", 10)))

# Check to make sure everything looks correct: most manuall assigned
model.multi.plot

ggplot(data = model.multi.plot, aes(x = Predictor, y = estimate, color = Model)) +
	geom_hline(yintercept = 0, lwd = 1.25, color = "grey80") + 
	geom_point(size = 3, position = position_dodge(width = 0.5)) +
	scale_color_manual(values = c("#C994C7", "#DF65B0", "#67001F"), 
		breaks = c("Neither -> Both", "Neither -> Just Crypto", "Neither -> Just Stocks")) +
	geom_errorbar(aes(ymin = lower.95, ymax = upper.95), lwd = 0.5, width = 0, position = position_dodge(width = 0.5)) +
	geom_errorbar(aes(ymin = lower.90, ymax = upper.90), lwd = 1.1, width = 0, position = position_dodge(width = 0.5)) +
	facet_wrap(~ Category, scales = "free_y") +
	labs(y = "Estimate (All Predictors Scaled from 0-1)") +
	coord_flip() +
	theme_bw()
dev.off()







###############################################################################################
# Table 1. Sample demographics and personality of full sample and currency owners
###############################################################################################

# educ == 4: associates
data$college.degree.plus <- ifelse(data$educ > 4, 1, 0)
table(data$college.degree.plus, data$educ)

# dummy top and bottom categories of religion importance for demographics
data$relig.not.import <- ifelse(data$religimportance == 1, 1, 0)
data$relig.very.import <- ifelse(data$religimportance == 4, 1, 0)
table(data$relig.not.import, data$religimportance)
table(data$relig.very.import, data$religimportance)

# dummy top and bottom categories of population density for demographics
data$rural <- ifelse(data$popdensity == 1, 1, 0)
data$urban <- ifelse(data$popdensity == 5, 1, 0)
table(data$rural, data$popdensity)
table(data$urban, data$popdensity)

## Demographic characteristics
# Full sample
demo.overall.sample <- data %>%
	select(age, male, white, black, hispanic, asian, college.degree.plus, married, 
		relig.very.import, relig.not.import, income, 
		stockowner, rural, urban, bornagain, catholic) %>%
	summarise_all(mean, na.rm = TRUE) %>%
	t()

# Cryptoowners
demo.crypto.sample <- data %>%
	select(age, male, white, black, hispanic, asian, college.degree.plus, married, 
		relig.very.import, relig.not.import, income, 
		stockowner, rural, urban, bornagain, catholic, cryptoowner) %>%
	filter(cryptoowner == 1) %>%
	select(!cryptoowner) %>%
	summarise_all(mean, na.rm = TRUE) %>%
	t()

# Cryptoowners
demo.diff.means <- data %>%
	select(age, male, white, black, hispanic, asian, college.degree.plus, married, 
		relig.very.import, relig.not.import, income, 
		stockowner, rural, urban, bornagain, catholic) %>%
	map_df(~ tidy(t.test(. ~ data$cryptoowner)), .id = 'var')


## Personality characteristics
# Full sample
personality.overall.sample <- data %>%
	select(openexp, conscientiousness, extraversion, agreeableness, emostable, 
		narcissism, machiavellianism, psychopathy, n2e, n4c) %>%
	summarise_all(mean, na.rm = TRUE) %>%
	t()

# Cryptoowners
personality.crypto.sample <- data %>%
	select(openexp, conscientiousness, extraversion, agreeableness, emostable, 
		narcissism, machiavellianism, psychopathy, n2e, n4c, cryptoowner) %>%
	filter(cryptoowner == 1) %>%
	select(!cryptoowner) %>%
	summarise_all(mean, na.rm = TRUE) %>%
	t()

# Cryptoowners
personality.diff.means <- data %>%
	select(openexp, conscientiousness, extraversion, agreeableness, emostable, 
		narcissism, machiavellianism, psychopathy, n2e, n4c) %>%
	map_df(~ tidy(t.test(. ~ data$cryptoowner)), .id = 'var')


rbind(
cbind(round(demo.overall.sample, digits = 4), round(demo.crypto.sample, digits = 4), round(demo.diff.means$p.value, digits = 3)),
cbind(round(personality.overall.sample, digits = 4), round(personality.crypto.sample, digits = 4), round(personality.diff.means$p.value, digits = 3))
)




###############################################################################################
# Table 2. Sample politics of full sample and currency owners
###############################################################################################

## Demographic characteristics
# Full sample
politics.overall.sample <- data %>%
	select(spendserve, trustgov, conspiratorialthinking, personalinflation, inflationexpectations, 
		libertarianvalues, pid7sr, symbideologycons) %>%
	summarise_all(mean, na.rm = TRUE) %>%
	t()

# Cryptoowners
politics.crypto.sample <- data %>%
	select(spendserve, trustgov, conspiratorialthinking, personalinflation, inflationexpectations, 
		libertarianvalues, pid7sr, symbideologycons, cryptoowner) %>%
	filter(cryptoowner == 1) %>%
	select(!cryptoowner) %>%
	summarise_all(mean, na.rm = TRUE) %>%
	t()

# Cryptoowners
politics.diff.means <- data %>%
	select(spendserve, trustgov, conspiratorialthinking, personalinflation, inflationexpectations, 
		libertarianvalues, pid7sr, symbideologycons) %>%
	map_df(~ tidy(t.test(. ~ data$cryptoowner)), .id = 'var')

# 
cbind(round(politics.overall.sample, digits = 4), round(politics.crypto.sample, digits = 4), round(politics.diff.means$p.value, digits = 3))



###############################################################################################
# Table A.1. Tabular results of base and full model from Figure 2.
###############################################################################################

the.table <- NA

for(i in 1:length(the.full.predictors)) {
	if(the.full.predictors[i] %in% base.model.plot$Predictor) {
		base.est <- round(base.model.plot$estimate[base.model.plot$Predictor == the.full.predictors[i]], digits = 3)
		base.se <- paste0("(", round(base.model.plot$std.error[base.model.plot$Predictor == the.full.predictors[i]], digits = 3), ")")
		base <- paste0(base.est, " ", base.se)
	} else {
		base <- "\\multicolumn{2}{c}{-}"
	}
	full.est <- round(full.model.plot$estimate[full.model.plot$Predictor ==  the.full.predictors[i]], digits = 3)
	full.se <- paste0("(", round(full.model.plot$std.error[full.model.plot$Predictor ==  the.full.predictors[i]], digits = 3), ")")
	the.row <- paste0(paste(the.full.predictors[i]), "&", base, "&", full.est, " ", full.se, "\\")
	the.table <- rbind(as.matrix(the.table), as.matrix(the.row))
}

the.table

###############################################################################################
# Table A.2. Tabular results of multinomial logit from Figure 4.
###############################################################################################
apply(
	cbind(as.matrix(paste0(model.multinom.predictors, " & ")), 
		as.matrix(paste0(round(coef(model.multinom)[1,], digits = 3), " (", round(sqrt(diag(vcov(model.multinom)))[1:29], digits = 3), ") & ")), 
		as.matrix(paste0(round(coef(model.multinom)[2,], digits = 3), " (", round(sqrt(diag(vcov(model.multinom)))[30:58], digits = 3), ") & ")), 
		as.matrix(paste0(round(coef(model.multinom)[3,], digits = 3), " (", round(sqrt(diag(vcov(model.multinom)))[59:87], digits = 3), ") \\")) 
	), 
	2, paste0
)


###############################################################################################
# Table A.3. Tabular results of regressions by PID
###############################################################################################

dem.model <- glm(cryptoowner ~ age + male + black + hispanic + asian + educ + married + popdensity +
								bornagain + catholic + religimportance + 
								income + stockowner + personalinflation + inflationexpectations +
								spendserve + trustgov + conspiratorialthinking + symbideologycons + libertarianvalues +
								openexp + conscientiousness + extraversion + agreeableness + emostable + 
								n2e + n4c + narcissism + machiavellianism + psychopathy, 
								family = binomial(link = "logit"), data = data01, subset = pid7sr < 0.5)

indp.model <- glm(cryptoowner ~ age + male + black + hispanic + asian + educ + married + popdensity +
								bornagain + catholic + religimportance + 
								income + stockowner + personalinflation + inflationexpectations +
								spendserve + trustgov + conspiratorialthinking + symbideologycons + libertarianvalues +
								openexp + conscientiousness + extraversion + agreeableness + emostable + 
								n2e + n4c + narcissism + machiavellianism + psychopathy, 
								family = binomial(link = "logit"), data = data01, subset = pid7sr == 0.5)

rep.model <- glm(cryptoowner ~ age + male + black + hispanic + asian + educ + married + popdensity +
								bornagain + catholic + religimportance + 
								income + stockowner + personalinflation + inflationexpectations +
								spendserve + trustgov + conspiratorialthinking + symbideologycons + libertarianvalues +
								openexp + conscientiousness + extraversion + agreeableness + emostable + 
								n2e + n4c + narcissism + machiavellianism + psychopathy, 
								family = binomial(link = "logit"), data = data01, subset = pid7sr > 0.5)

cbind(round(dem.model$coefficients, digits = 3), paste0("(", round(sqrt(diag(vcov(dem.model))), digits = 3), ")"),
	round(indp.model $coefficients, digits = 3), paste0("(", round(sqrt(diag(vcov(indp.model))), digits = 3), ")"),
	round(rep.model $coefficients, digits = 3), paste0("(", round(sqrt(diag(vcov(rep.model))), digits = 3), ")")
)















